Contents
function hw4() clear(); close all; load('velb4degl.mat'); % Figure out how to partition this to get 50 DOF dof = 50; len_complete = size(velb4, 2); points_per = floor(len_complete / dof); for i=1:dof records.velocity(:,:,i) = velb4(:,(i-1)*points_per+1:i*points_per); records.time = 0:0.625:0.625*points_per; % times records.space = 0:2:2*219; % spacing end % Task 1: Present power spectral estimates of the horizontal velocity % Ehat(sigma) at 50 degrees % of freedom for two ranges separated by ?r=24m. Plot log-log on the same plot. Identify % the swell peak, wind-wave peak. State the frequency in cycles/s, please. % Get two records, 24m apart. % From previous problems, we know that data in the spatial range 10:201 is % good. Each row is separated by 2m. So, pick two rows 12 indexes apart. % Index 10 = 20m (or maybe more... don't know what index 1 is, actually.) % index 22 = 44m %velocity_close(:,:) = records.velocity(10, :, :); %velocity_far(:,:) = records.velocity(22, :, :); % Get two records, 24m apart. % From previous problems, we know that data in the spatial range 10:201 is % good. Each row is separated by 2m. So, pick two rows 12 indexes apart. % Index 100 = 200m (or maybe more... don't know what index 1 is, actually.) % index 112 = 124m % Both records will use the same frequency spacing, so we don't need % different variables for that. sample_interval = 0.625; % seconds velocity_close = reshape(velb4(100,:), [], 1); velocity_far = reshape(velb4(112,:), [], 1); [ehat_close_real, ehat_close, freq_bins_real, freq_bins, delta_f] = ... varspec_est_multidof(velocity_close, sample_interval, dof); [ehat_far_real, ehat_far, freq_bins_real, freq_bins, delta_f] = ... varspec_est_multidof(velocity_far, sample_interval, dof); figure(); clf(); loglog(freq_bins_real, ehat_close_real); hold on; loglog(freq_bins_real, ehat_far_real, 'r--'); % I don't think we have a single axis label here. %ylabel('Variance of velocity (m^2/s * hours)') xlabel('Frequency (cycles/sec)'); legend('Close Range (200m)','Far Range (224m)'); title('Spectral density estimates using two ranges'); ylabel('Variance spectral density of velocity (cm/s) ^-^2'); % Task 2: Present co and quadrature spectral estimates at 50 d.o.f. Plot log (frequency, % cps) vs. linear ordinate. [cross_spec, freq_bins] = crossspectrum(velocity_close, velocity_far, sample_interval, dof); figure() % I co = real(cross_spec(1:60)); semilogx(freq_bins, co); hold on; % Q quad = imag(cross_spec(1:60)) semilogx(freq_bins, quad, 'r--'); % Task 3: Present estimates of coherence and phase (in degrees) at 50 d.o.f. % On the phase plot, re-plot the phase estimate in each frequency band as an asterisk (*), % but only in those frequency bands where the coherence is significant at the 67% % confidence level. coherence = abs(cross_spec(1:60)) ./ sqrt(ehat_close_real .* ehat_far_real); % angle(z) is the same as atan2(imag(z), real(z)) phase_degrees = angle(cross_spec(1:60)) * (180/pi); figure() semilogx(freq_bins, abs(cross_spec(1:60))); figure() semilogx(freq_bins, phase_degrees); end function [cross_spec, freq_bins] = crossspectrum(data_a, data_b, sample_interval, dof) % We get 2 DOF from each set of data (thanks to real/imag components. num_sets = ceil(dof / 2); len_complete = size(data_a, 1); points_per = floor(len_complete / num_sets); delta_f = (1/sample_interval) / points_per; % same as fs / N freq_bins = 0:delta_f:(points_per/2 * delta_f); % Get the variance spectral densitity for each subset. for i=1:num_sets start_idx = (i-1) * points_per + 1; end_idx = start_idx + points_per - 1; % normalize FFT by n^2 a_a(:,i) = fft(data_a(start_idx:end_idx)) ./ points_per^2; a_b(:,i) = fft(data_b(start_idx:end_idx)) ./ points_per^2; crosses(i,:) = bsxfun(@times, conj(a_a(:, i)), a_b(:,i)); end % Get the average values (average of columns) cross_spec = mean(crosses) ./ delta_f; % get rid of the first (mean) value. cross_spec = cross_spec(2:end); freq_bins = freq_bins(2:end); end
Variance Spectral Density Functions
Take the given data and provide a variance density estimate The freq_bins array will be in reciprocal units to the units of sample_interval. ehat has units of (data_units^2) per (1/sample_interval_units)
function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ... varspec_est_real(data, sample_interval) % First, do the FFT to get an array of Fourier coefficients a = fft(data); % Figure out our frequency bins. N = length(data); delta_f = (1/sample_interval) / N; % same as fs / N freq_bins_real = 0:delta_f:(N/2 * delta_f); freq_bins = 0:delta_f:(N * delta_f); % We are dealing with a real signal, so we only care about positive, real % freqency components. To get the power, we want to "fold" both sides of % the spectrum together. Since they are the same, we can just take one % side and multiply by 2. % To preserve variance, we have to normalize this result by dividing by % N^2, due to the Matlab implementation of fft. ehat_real = 2 * abs(a(1:length(freq_bins_real))).^2 ./ delta_f ./ ... (N^2); % Also, figure out the full spectrum ehat = abs(a(1:length(freq_bins)-1)).^2 ./ delta_f ./ (N^2); % The first frequency bin will hold the mean, not a component of the % variance. Since we are interested in the variance spectrum, drop % the first bin. ehat = ehat(2:end); ehat_real = ehat_real(2:end); freq_bins = freq_bins(2:end); freq_bins_real = freq_bins_real(2:end); end function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ... varspec_est_avg(data, sample_interval, set_size, number_of_sets) % Get the variance spectral densitity for each subset. for i=1:number_of_sets start_idx = (i-1) * set_size + 1; end_idx = start_idx + set_size - 1; [ehat_parts_real(i, :), ehat_parts(i, :), freq_bins_real, freq_bins, delta_f] = ... varspec_est_real(data(start_idx:end_idx), sample_interval); end % form the average % Get the average values (average of columns) ehat = mean(ehat_parts); ehat_real = mean(ehat_parts_real); end function [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ... varspec_est_multidof(data, sample_interval, dof) % We get 2 DOF from each set of data (thanks to real/imag components. num_sets = ceil(dof / 2); len_complete = size(data, 1); points_per = floor(len_complete / num_sets); [ehat_real, ehat, freq_bins_real, freq_bins, delta_f] = ... varspec_est_avg(data, sample_interval, points_per, num_sets); end
quad = Columns 1 through 16 -0.0017 0.0011 -0.0039 0.0021 0.0584 0.0736 0.0226 0.0163 0.0142 0.0103 0.0122 0.0425 0.0374 0.0372 0.0061 0.0019 Columns 17 through 32 -0.0021 -0.0049 -0.0024 -0.0050 -0.0068 -0.0060 -0.0011 0.0003 0.0037 0.0023 0.0016 0.0012 0.0022 0.0012 0.0002 0.0018 Columns 33 through 48 -0.0003 -0.0006 0.0007 0.0003 -0.0009 0.0007 -0.0012 -0.0003 -0.0002 0.0006 0.0006 0.0004 0.0002 0.0001 -0.0010 -0.0010 Columns 49 through 60 -0.0010 0.0023 -0.0012 -0.0020 0.0006 0.0008 -0.0009 0.0008 0.0007 -0.0004 -0.0009 0



